function FilDyn = CKFILTER (Dyn, Span, Kfb, W, PI)

%--------------------------------------------------------------------------
% function CKFILTER *Created by PEIWEN CONG* 23.11.10
% CK forward backward non-linear filter implementation
% reference doi:10.1016/0165-0270(91)90118-J
%
%##########NOTE:Boundary Conditions###########
% See CKMEANARRAY for details
% Set to 0 if no points
% Set to mean(short) if not long enough
%#############################################
%
% INPUT
% Dyn:      time series of observations
% Span:     analysis window size
% Kfb:        forward (backward) predictor banks
% PI:       prior estimated probabilities
% W:        weight factors
%
% OUTPUT
% FilDyn:   filtered dynamics
%
% FUNCTIONS LOADED
% CKMEANARRAY
%--------------------------------------------------------------------------

% input check
nbank = length(Kfb);
Kf = Kfb(1:nbank/2);
Kb = Kfb(nbank/2+1:nbank);

if nargin < 5
    PI = 1/nbank;
end

% obtain K forward-backward vector: Kfb
% e.g. [ 16 8 4 2 -2 -4 -8 -16 ]

% obtain weight factors vector: W
if length(W) ==1 
    W = W.*ones(1,nbank);
end

%disp(W)

% obtain prior estimation vector: PI
if length(PI) ==1 
    PI = PI.*ones(1,nbank);
end

%disp(PI)

% initialization
FilDyn = zeros(1,length(Dyn));

% loop scan Dyn
for k=1:length(Dyn)
    KfSpan = [k-Span+1:k];
    for i = 1:Span
        if (KfSpan(i) < 1)
            KfSpan(i) = 1;
        end
    end
    KbSpan = [k:k+Span-1];
    for i = 1:Span
        if (KbSpan(i) > length(Dyn))
            KbSpan(i) = length(Dyn);
        end
    end
    
    %SDArray (M * 2K)
    SDfArray = kron(Dyn(KfSpan),ones(nbank/2,1));
    SDbArray = kron(Dyn(KbSpan),ones(nbank/2,1));
    SDArray = [SDfArray; SDbArray];
    %FBArray (M * 2K)
    FBfArray = CKMEANARRAY(Dyn,KfSpan,Kf);
    FBbArray = CKMEANARRAY(Dyn,KbSpan,Kb);
    FBArray = [FBfArray; FBbArray];
    %disp(sum((SDArray-FBArray).^2,2)'.^(-W))
	weight = PI.*(sum((SDArray-FBArray).^2,2)'.^(-W));
    % Inf weight for boundary conditions
    if sum(isinf(weight)) ~= 0
        weight = isinf(weight);
    end
    normalweight = weight ./ sum(weight);
    %disp(normalweight)
    %single k
	FBArrayk = CKMEANARRAY(Dyn,k,Kfb);
    %disp(FBArrayk)
    FilDyn(k) = normalweight*FBArrayk;
end

